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20.  Continued 

significantly  larger  time  steps  with  comparable  accuracy.  In  many,  if  not  most,  problems  it  is 
not  obvious  which  integration  method  is  more  efficient.  In  this  study  the  Ncwmark  Beta- 
Method  is  examined  for  stability,  accuracy,  and  efficiency,  wherein  Beta  = 0 provides  an 
explicit  algorithm,  while  Beta  SO  provides  an  implicit  algorithm.  Both  algorithms  arc  used  in 
the  same  finite  element  program  to  solve  a soil-structure  boundary  value  problem  composed 
of  a cylindrical  steel  shell  encased  in  a relatively  soft  rock-like  material  and  subjected  to  a 
surface  blast  loading.  Bor  this  problem  with  linear  system  properties,  the  implicit  method 
was  significantly  more  efficient  as  measured  by  computer  time,  bor  nonlinear  systems,  the 
two  methods  arc  approximately  equivalent  in  efficiency.  A combined  explicit-implicit 
integration  technique  is  proposed  for  these  types  of  interaction  problems  with  two  or  more 
materials.  The  combined  explicit-implicit  algorithm  employs  explicit  integration  in  the  soft 
material  and  implicit  integration  in  the  stiff  material  with  a potential  increase  in  efficiency 
by  an  order  of  magnitude  over  cither  method  applied  individually. 
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Direct  integration  techniques  (step-by-step)  arc  widely  used  for  the  time  integration  of  dis- 
cretized equations  of  motion  that  result  from  applying  numerical  techniques  such  as  the  finite  ele- 
ment method  to  structural  dynamic  problems.  In  this  study,  the  Ncwmark  Beta-Method  is 
examined  for  stability,  accuracy,  and  efficiency,  wherein  Beta  a o provides  an  explicit  algorithm,/  • 
while  Beta  * 0 provides  an  implicit  algorithm.  Both  algorithms  arc  used  in  the  same  finite  element 
program  to  solve  a soil-structure  boundary  value  problem  composed  of  a cylindrical  steel  sbccl  en- 
cased in  a relatively  soft  rock-like  material  anil  subjected  to  a surface  blast  loading.  The  implicit 
method  was  significantly  more  efficient  as  measured  by  computer  time.  For  nonlinear  systems,  the 
two  methods  were  approximately  equivalent  in  efficiency.  A combined  explicit-implicit 
integration  technique  is  proposed  that  employs  explicit  integration  in  the  soft  material  and 
implicit  integration  in  the  stiff  material  with  a potential  increase  in  efficiency  by  an  order  of 
magnitude  over  either  method  applied  individually. 
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INTRODUCTION 

In  the  past  two  decades  the  finite  element  method  has  been  widely 
used  for  solving  structural  dynamics  problems.  The  numerical  approxima- 
tions inherent  in  the  finite  element  technique  include  two  broad  cate- 
gories: (a)  spatial  discretization  and  (2)  temporal  discretization.  The 

former  category  deals  with  the  approximation  errors  and  efficiency  of  the 
assumed  spatial  distribution  of  the  primary  variables  (e.g. ,,  displacement 
field),  while  the  latter  category  deals  with  the  approximation  errors  and 
efficiency  involved  with  the  time  integration  of  the  discretized  equations 
of  motion.  In  this  study  only  the  latter  category  is  pursued  and  the 
discussion  is  restricted  to  direct  integration  methods  (i.e. , step-by- 
step).  Other  methods  of  solving  the  equations  of  motion  such  as  eigen- 
value analysis  or  method  of  characteristics  are  not  considered  because  of 
their  limitations  for  nonlinear  problems. 

Direct  integration  schemes  are  generally  classified  in  two  groups: 
implicit  methods  and  explicit  methods.  This  report  examines  these  two 
methods  in  the  light  of  numerical  stability,  accuracy,  and  efficiency  for 
a class  of  problems,*  characterized  by  a cylindrical  structural  liner 
embedded  in  a homogeneous  rock-like  material  and  subjected  to  a surface 
blast  loading  [1,  2], 

Background 

Numerous  implicit  and  explicit  integration  operators  are  in  current 
use  and  have  been  reported  extensively  in  engineering  literature.  Some 
of  the  more  popular  implicit  integration  schemes  are:  the  Newmark 

B-method  (3  h 0)  [3],  the  Houbolt  method  [4],  and  the  Wilson  0-method  [5], 
Explicit  schemes  are  more  commonly  referred  to  by  a finite  difference 
operator  such  as  the  second  central  difference  scheme  [6],  Note  the 
Newmark  B^method  with  S ■ 0 is  equivalent  to  the  second  central  difference 
3cherne  and  is  an  explicit  method. 

The  fundamental  difference  between  an  explicit  and  implicit  integra- 
tion scheme  is  that  in  the  explicit  case  the  displacement  vector  for  the 
current  time  step  can  be  predicted  directly  from  the  known  displacements, 
velocities,  and  accelerations  of  the  previous  time  steps,  whereas  in  the 
implicit  case  the  current  displacements  are  related  to  current  accelerations 
(and  possibly  current  velocities)  as  well  as  responses  from  previous  time 
steps.  As  a result,  the  solution  algorithm  for  the  implicit  scheme  requires 
assembling  a global  mass  and  stiffness  matrix  into  a set  of  coupled  alge- 
braic equations  and  solving  the  system  by  some  technique  such  as  Gaussian 
elimination. 


*0f  current  interest  to  the  Defense  Nuclear  Agency  (DNA). 
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On  the  other  hand,  the  explicit  algorithm  only  requires  assembling 
the  mass  matrix  and  solving  for  the  current  acceleration.  The  restoring 
force  vector  (stiffness  matrix  times  predicted  displacement  vector)  is 
formed  on  the  right-hand  side  at  the  element  level,  thereby  eliminating 
the  need  of  assembling  the  global  stiffness  matrix.  If  a lumped  mass 
procedure  is  used  (as  is  generally  done) , the  accelerations  are  uncoupled 
and  can  be  computed  rapidly. 

The  trade-off  one  pays  for  the  rapid  .explicit  algorithm  is  the 
requirement  of  a small  time  step  to  avoid  numerical  instability.  When 
instability  occurs  the  computed  responses  become  wildly  erratic,  oscil- 
lating many  orders  of  magnitude.  Numerical  instability  is  governed  by 
the  frequency  content  of  the  discretized  system  and  will  be  discussed 
later. 

The  potential  advantage  of  some  implicit  integration  schemes  is  that 
they  are  unconditionally  stable,  allowing  a much  larger  time  step  con- 
trolled by  integration  accuracy  (as  opposed  to  integration  stability) . 

In  short,  if  computational  efficiency  is  measured  by  the  computer 
cost  (time)  to  accurately  integrate  the  equations  of  motion  over  a given 
time  range  of  interest,  then  the  explicit  method  has  the  advantage  of 
rapid  solutions  per  time  step  but  the  disadvantage  of  requiring  many  time 
steps.  The  reverse  is  true  for  the  implicit  method. 

It  cannot  be  categorically  claimed  that  one  method  is  more  efficient 
than  the  other.  Rather,  the  question  of  efficiency  is  highly  problem 
dependent,  involving  such  factors  as:  finite  element  mesh  topology, 

types  of  elements,  loading  function,  material  properties,  degree  of  non- 
linearity, method  of  characterizing  mass  matrix,  and  the  size  of  the 
system.  Discussion  of  these  influences  are  given  in  References  7 and  8. 

Recently  the  idea  has  been  advanced  to  simultaneously  employ 
implicit  and  explicit  integration  operators  over  different  regions  of 
the  finite  element  mesh  to  take  optimum  advantage  of  each  method.  The 
above  notion,  as  well  as  explicit  versus  implicit  efficiency  studies,  is 
pursued  in  the  course  of  this  investigation. 

Objective  and  Scope 

This  investigation  is  restricted  to  the  family  of  integration 
operators  contained  in  Newmark’s  8-method  [3]  wherein  8 * 0 provides  an 
explicit  scheme  and  8 **  0 provides  an  implicit  scheme.  The  test  boundary 
value  problem  to  be  investigated  is  an  elastic  steel  cylinder  encased  in 
a homogeneous,  elastic,  rcck-like  material  subjected  to  a surface  blast 
loading. 

Within  the  above  framework,  the  objectives  of  this  study  include: 

(1)  Identification  and  discussion, of  explicit  and  implicit  algorithms, 
including  considerations  of  numerical  stability  and  accuracy. 

(2)  Comparison  of  efficiency  between  explicit  and  implicit  for  the 
test  boundary  value  problem  and  estimates  for  a nonlinear  case. 

(3)  Development  of  a combined  explicit-implicit  integration 
algorithm  for  superior  efficiency. 
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Approach 


It  is  recognized  that  an  efficiency  comparison  of  explicit  versus 
implicit  schemes  (Objective  2)  is  relative  and  subjective.  That  is,  the 
relative  merits  of  the  two  integration  schemes  are  not  only  dependent  on 
the  nature  of  the  boundary  value  problem,  but  also  on  the  programming 
skill  of  the  program  authors  and  the  particular  costing  algorithm  of  the 
host  computer  site. 

To  minimize  the  influence  of  different  programmers,  the  finite 
element  program  FEAP  [9]  was  adopted  for  this  study.  FEAP  js  the  only 
known  general  purpose  finite  element  code  containing  both  explicit  and 
implicit  algorithms  based  on  the  Newmark  6-method.  It  is  a linear  code 
(also  linear  viscoelastic)  with  two-  and  three-dimensional  elements. 

With  regard  to  the  host  computer  site  (CDC  6600) , only  central  pro- 
cessing time  will  be  considered  for  efficiency  comparisons  without  regard 
to  input/output  (I/O)  time.  Altnough  I/O  time  may  be  significant,  partic- 
ularly for  out-of-core  solution  methods,  it  is  not  considered  in  this 
study  because  the  associated  costs  are  arbitrary  and  dependent  on  the 
costing  algorithm  of  the  operating  site. 

In  the  following  pages,  the  Newmark  6-method  is  reviewed  in  the 
context  of  explicit  and  implicit  algorithms,  together  with  a discussion 
of  numerical  stability  and  integration  accuracy.  Next,  the  FEAP  code  is 
used  to  ascertain  the  efficiency  of  explicit  versus  implicit  for  the 
stated  test  problem.  Lastly,  a combined  explicit-implicit  algorithm  is 
presented  followed  by  conclusions  and  recommendations. 


THEORETICAL  CONSIDERATIONS  OF  NEWMARK  6-METHOD 


For  this  investigation,  we  consider  the  spatially  discretized  matrix 
equation  of  structural  dynamics  given  by: 

Mu  + Ku  - f (1) 


where : 


M - mass  matrix 
K « stiffness  matrix 

3B 

u * acceleration  vector 
u « velocity  vector 
u ■ displacement  vector 
f ■ external  load  vector 
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Equation  1 is  a coupled  system  of  second-order  ordinary  differential 
equations  which  represents  an  initial-value  problem  for  finding  the  dis- 
placement ^unction  u(t)  satisfying  Equation  1 and  also  satisfying  the 
initial  conditions  u(0)  = uD  and  u(0)  * u0,  where  u0  and  u0  are  given 
initial  data.  The  matrix  is  assumed  constant,  symmetric,  and  positive 
definite.  The  stiffness  matrix  K is  also  assumed  positive  definite,  but 
not  necessarily  constant,  to  reflect  a changing  stiffness  due  to  material 
nonlinearity. 

Newmark  6-Method 

A step-by-step  operator  suitable  for  numerically  integrating 
Equation  1 is  given  by  Newmark  [3]  as: 

Ht«t  ’ St  + itt  (1  - y)ut  + Yut+it]  (2) 

-t-Mt  " -t  + 4t«t  + *t2I  C /2  " 6>ut  + 6ut+4t)  (3) 


where:  At  - time  step 

6 - Newmark  B-parameter,  0 5 B - 1/4 
y * Newmark  numerical  damping  parameter  (y  * 1/2) 

The  damping  parameter  y introduces  numerical  damping  when  y * 1/2.  In 
this  study,  y is  taken  as  1/2  to  avoid  complications  of  numerical  damping. 
Subscripts  denote  the  time  at  which  the  response  is  evaluated. 

Specification  of  the  B-parameter  (B  * 0)  leads  to  a variety  of  well- 
known  implicit  schemes.  For  example,  B * 1/4  is  equivalent  to  the  average 
acceleration  scheme,  and  B **  1/6  is  equivalent  to  the  linear  acceleration 
scheme.  For  the  special  case  6 ■ 0,  an  explicit  scheme  results,  as  can 
be  observed  from  Equation  3 wherein  displacements  at  time  t + At  are  pre- 
dicted from  known  responses  at  the  previous  time  t. 

Computational  algorithms  for  implicit  (B  *»  0)  and  explicit  (B  * 0) 
integrations  are  outlined  in  the  next  sections. 

Implicit  Algorithm 

It  is  assumed  the  responses  ut,  ut,  and  ut  are  known  for  time  t,  and 
the  objective  is  to  determine  these  responses  at  time  t + At.  To  this 
end,  the  following  three  steps  constitute  an  implicit  algorithm: 
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Step  1.  From  Equation  3,  express  Ut+At  as  a function  of  ut+^t  and 
previous  responses  to  get: 


St+it  " <,/64t  >2t+At  • st 


where : afc  * 


Step  2.  Inserting  the  above  expression  for  Ut+At  into  Equation  1, 

ut+^t  is  determined  from  the  solution  of  the  coupled  system: 


iSAt 


2)  — + * 


*t+At 


-t+At  + -~t 


Step  3.  Update  accelerations  using  the  expression  in  Step  1,  and 
update  velocities  using  Equation  2,  i.e. , 


~t+At 


1 


,8At 


2 1 -t+At  ' -t 


At 


St+At  ■ St  + T (Bt  + “t+At1 


The  above  algorithm  is  repeated  successively  for  each  time  step  through 
the  time  range  of  interest. 

Step  2 of  the  above  algorithm  requires  assembling  the  global  mass 
and  stiffness  matrix  and  solving  the  coupled  system  by  techniques  such 
as  elimination  (i.e.,  triangularization) , which  is  a costly  operation. 

If  the  system  is  linear  and  a constant  size  time  step  is  used,  the  tri- 
angularization can  be  done  once  and  for  all  at  the  outset  of  the  proce- 
dure so  that  subsequent  steps  only  require  modifying  the  right-hand  side 
and  performing  back  substitution.  However,  for  nonlinear  systems  the 
stiffness  matrix  may  have  to  be  reformed  and,  hence,  resolved  at  every 
time  step  and  perhaps  iterations  within  the  time  step,  depending  on  the 
degree  of  nonlinearity  and  the  methodology  employed.  This  disadvantage 
of  the  implicit  method  is  not  shared  by  the  explicit  algorithm  discussed 
next. 


Explicit  Algorithm 

As  before,  it  is  assumed  the  responses  ut,  ut>  and  iit  are  known  and 
the  objective  is  to  determine  the  responses  at  time  t+At.  The  following 
three  steps  constitute  an  explicit  algorithm. 
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Step  1.  From  Equation  3 with  6*0,  the  displacements  uL+£t  are 

predicted  directly  as  a function  of  the  previous  responses: 


-t+At 


*t 


where:  dt 


2 

. . .At  .. 
u.  + Atu„  + — — u 

— t “t  X.  — t 


Step  2.  Inserting  the  above  expression  for  ut+^t  into  Equation  1 , 
the  accelerations  ut+^t  are  given  by  the  solution  of: 


= -t+At  “ -t+At  “ = ~t 


Step  3.  Lastly,  the  velocities  are  updated  from  Equation  2, 

i.e. : 


• • i At  r"  i ••  i 

-t+At  ‘ St  + T tit  + -t+At 


Unlike  Step  2 of  the  implicit  algorithm,  Step  2 of  the  explicit 
algorithm  only  requires  assembling  the  global  mass  matrix  If  the 
lumped  mass  technique  is  adopted,  then  is  diagonal  and  the  solution 
for  ut+At  is  trivial  (otherwise  £f  must  be  triangularized  once  at  the 
outset).  Furthermore,  the  internal  force  vector  represented  by  Kdt  can 
be  easily  computed  at  the  element  level  circumventing  the  need  of  formally 
assembling  the  global  stiffness  K.  Accordingly,,  there  is  no  significant 
penalty  in  the  algorithm  for  nonlinear  systems  since  changes  in  & are 
easily  accommodated  by  right-hand-side  operations  at  the  element  level. 

The  apparent  advantages  of  the  explicit  over  implicit  algorithm  are 
mitigated  by  the  concern  of  numerical  stability  which  limits  the  size  of 
the  integration  step  At.  These  notions  are  discussed  next. 

Stability  of  Newmark  6-Method 

Numerical  stability  of  a solution  algorithm  may  generally  be  assured 
by  proper  choice  of  the  integration  step  size.  Since  stability  increases 
with  decreasing  step  size,  a sufficiently  small  step  can  generally  be 
found  to  provide  a stable  solution.  Therefore,  the  question  to  be  explored 
is,  what  size  time  step  is  required  for  stability  of  the  Newmark  6-method. 
The  following  development  is  restricted  to  linear  systems.  Stability  for 
nonlinear  systems  has  been  addressed  on  an  energy  basis  [IQ,  11]. 

It  is  well-known  that  a linear  system  represented  by  Equation  1 can 
be  transformed  into  an  equivalent  set  of  uncoupled  equations  through 
modal  analysis.  A typical  uncoupled  modal  equation  may  bp  written  as: 


6 


X 


a)2  X 


P(t) 


(4) 


+ 


where  X is  the  transformed  displacement  called  the  normal  coordinate, 
u is  the  associated  frequency  of  vibration,  and  p(t)  is  the  transformed 
loading  function. 

It  is  asserted  that  the  numerical  integration  of  all  the  modal 
equations  (typified  by  Equation  4)  using  the  same  At  is  equivalent  to 
the  numerical  integration  of  the  coupled  system  (Equation  1).  Accordingly, 
numerical  stability  can  be  simplified  to  the  investigation  of  Equation  4. 

To  this  end,  the  Newmnrk  integration  operator  (Equations  2 and  3) 
can  be  combined  and  expressed  independently  of  velocities  by  taking  the 
difference  of  Equation  3 at  time  t + At  from  its  value  at  time  t and 
replacing  the  resulting  velocity  difference  by  Equation  2.  This  gives: 


Xt+At  • 2Xt 


+ Xt-At  " At‘ ^Xt+At  + (1-26)Xt  + PXt-At] 


(5) 


T’here  the  nodal  displacement  u is  replaced  by  the  normal  coordinate  X. 

For  the  case  of  free  vibration  p(t)  * 0,  Equation  4 is  introduced 
into  the  integration  operator  (Equation  5)  at  times  t - At,  t,  and  t + At, 
resulting  in  the  difference  scheme: 


t+At 


bXt  + 


t-At 


(6) 


where 


- (mAh)2  (1  - 28) 
1 + (uAt)2  3 


(7) 


Equation  6 represents  a step-by-step  finite  difference  scheme  for 
successively  updating  the  displacements  Xt+^t  in  terms  of  the  previous 
displacements  Xt  and  Xt.^t, 

For  a given  set  of  initial  conditions,  the  free  vibration  of  a linear 
system  must  be  harmonic  and  bounded.  Thus,  the  stability  question  is 
stated:  for  a particular  $ and  maximum  frequency  ax,  how  large  a time 

step,  At,  may  be  taken  so  that  Equation  6 will  produce  bounded  3nd 
harmonic  results?  This  question  is  answered  by  determining  the  difference 
solution  which  is  achieved  by  initially  assuming  a solution  of  the  form: 

(t. /At) 

X.  - r (8) 

where  k is  the  time  step  counter  k * t^/At,  and  r is  to  be  determined. 
Inserting  Che  above  into  Equation  6 leads  to  the  characteristic  equation 
for  r: 
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- br 


+ 1 


whose  roots  r^ , ^ in  complex  polar  form  are: 


r1»  r2 


„±1{ 1 


(9' 


(10) 


where 


(ID 


iqtk/At 

Since  e * cr  s (qt^/At)  + i sin  (qt^/At) , the  difference 

solution  of  Equation  6 may  be  written  in  the  form: 


(12) 


where  A and  B are  real  constants  determined  from  initial  conditions; 

A « X(0),  B « (At/q)  X(0).  The  above  equation  yields  a bounded,  harmonic 
solution  for  Xtk,  providing  q is  real.  Clearly,  from  Equation  11,  q is 

real  if  (1  - b2/4)  J 0.  This  provides  the  stability  criterion,  such  that 
replacing  b with  Equation  7,  the  maximum  allowable  time  step  is  given  by: 


where  T is  the  shortest  period  of  vibration  of  the  system  given  by 
T * ^/“hnax* 

Evaluating  Equation  13,  the  allowable  time  step  for  the  explicit 
case,  B * 0,  is  At  S C.318  T.  For  implicit  examples,  B “ 1/6  gives 
At  S 0.551  T,  but  for  B ■ 1/4  there  is  no  finite  limit  on  At  for  stability. 
Hence,  8 ■ 1/4  provides  an  unconditionally  stable  operator. 

For  all  cases  B < 1/4,  the  allowable  time  step  to  insure  stability 
is  dependent  on  the  shortest  period  of  the  system  as  given  by  Equation  13. 
Unfortunately,  the  value  of  the  shortest  period  (or  highest  frequency)  is 
generally  not  known  and  troublesome  to  calculate.  As  an  alternative,  it 
is  often  more  convenient  to  determine  the  stability  limit  of  At  by  a 
heuristic  approach  given  next. 

Heuristic  Stability  Criterion 

In  lieu  of  using  Equation  13  for  selecting  a stable  time  step,  the 
following  relationship  may  be  used: 
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At*  < (h/c)min  >/ 1/(1  - 46) 


(14) 


where  h is  the  shortest  side  '.f  a finite  element,  c is  the  maximum  wave- 
speed  of  the  element  materia' , and  (h/c)min  implies  the  controlling 
element  in  the  mesh  where  t ■ ratio  h/c  is  a minimum.  (For  isotropic 
elastic  materials,  the  maximum  wavespeed  is  c * (E(1-v)/[p(1+v)  ( 1 -2v) ] }'1/2,  ' 
where  E * Young’s  Modul*  s,  v * Poisson’s  ratio,  and  p « mass  density.) 

The  genesis  of  Equation  14  is  based  on  heuristic  arguments  for  the 
explicit  algorithm  (8  * 0)  and  modified  by  the  factor  y 1 / ( 1 - 48)  for 
implicit  cases  as  discussed  in  the  following. 

Consider  the  continuum  body  shown  in  Figure  1 with  an  arbitrary 
finite  element  topology  drawn  on  the  body  and  focus  attention  on  the  node 
common  to  the  four  shaded  elements.  If  this  point  is  perturbed  by  an 
external  agency,  continuum  theory  requires  the  excitation  to  travel  with 
a sonic  wavespeed  c.  Therefore,  after  a time  interval,  At,  the  perturbed 
area  is  inscribed  by  a circle  of  radius  rc  - cAt,  shown  by  the  dashed 
line  in  Figure  1.  Now,  for  the  corresponding  explicit  algorithm,  the 
perturbed  area  for  one  time  step  can  be  no  greater  than  the  area  of 
numerical  coupling;  i.e.,  the  internal  forces  Kd  generated  by  an  excita- 
tion of  the  ■r.ommon  node  during  one  time  step  are  coupled  no  further  than 
the  nodes  of  the  shaded  elements.  This  condition  holds  for  each  and 
every  node  as  it,  in  turn,  becomes  excited.  Therefore,  in  order  to  insure 
that  the  “numerical  wavespeed’’  of  the  explicit  algorithm  can  excite  an 
area  at  least  as  large  as  the  area  of  sonic  travel,  it  is  required  that 
h i rc  or,  equivalently,  At  i (h/c)min. 

With  regard  to  the  implicit  algorithm  8 *»  0,  the  numerical  coupling 
is  greater  because  we  are,  in  effect,  operating  with  £"1  and  not  the 
banded  matrix  j£.  Since  K"  1 is  generally  fully  populated  (although 
weakly  coupled  between  distant  nodes),  the  area  of  numerical  coupling  is 
effectively  greater.  Accordingly,  the  maximum  allowable  time  step  is 
increased  by  the  factor  yi / ( 1 - 48) , which  is  the  ratio  from  Equation  13: 
At(B  \ 0)  v At(8  - 0), 

Experience  has  shown  that  Equation  14  provides  a good  estimate  of  an 
allowable  time  step  for  stability. 

Accuracy  of  Newmark  8-Method 

The  selection  of  At  to  insure  numerical  stability  is  a necessary 
condition  for  a meaningful  solution;  however,  it  is  not  automatically  a 
sufficient  condition  to  insure  that  the  numerical  results  are  a gooc 
approximation  of  the  original  differential  equation.  This  becomes  the 
question  of  accuracy.  An  obvious  example  is  the  case  8*1/4  which  has 
no  stability  limit  for  At  so  that  At  is  controlled  strictly  from  accuracy 
considerations. 

Accuracy  can  be  studied  by  comparing  the  exact  solution  of  the 
original  differential  equation  (Equation  4)  with  the  corresponding  dif- 
ference solution  given  by  Equation  12,  Specifically,  for  the  case  of 
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Figure  1.  Representation  of  numerical  and  continuum  wavespeed. 


free  vibration,  the  exact  solution  has  the  form  X * A cos  wt  + B sin  a>t, 
whereas  Equation.  12  has  the  form  - A cos  ut^  + B sin  wt^,  where 
5 ■ q/At.  For  harmonic  similarity,  the  numerical  frequency  ui  must  be  a 
good  approximation  of  the  actual  frequency  w;  or,  equivalently,  the 
numerical  period  T ■ 2r/u  must  be  a good  approximation  of  the  actual 
period  T - 2ir/w. 

The  ratio  of  the  approximate  period  to  the  actual  period  provides 
a measure  of  accuracy  [12].  With  the  aid  of  Equation  11,  this  ratio  may 
be  written  as: 


_T_  Zv  (At/T) 

T sin'*1  (/l  - b2/4) 

, 2 

2 “ ( 2rr  (1-2B) 

where  b 2 

1 + (2tt  ~)  3 

As  At  -*■  0,  T/T  1 , illustrating  that  the  approximate  solution 
approaches  the  exact  solution  in  the  limit.  However,  for  increasing 
values  of  At,  T/T  diverges  from  unity  on  a path  dependent  on  3.  Figure  2 
illustrates  this  trend  for  3*0,  1/6,  and  1/A  as  a function  of  At/T.  It 
is  observed  that  T/T  begins  to  significantly  deviate  from  unity  when  At/T 
is  larger  than  say  0.2,  which  is  smaller  than  the  stability  limits. 


05) 


(16) 
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From  the  above  observations,  it 
would  appear  that  accuracy  consider- 
ations place  a more  stringent  limit 
on  the  maximum  time  step  than  does 
the  stability  criterion.  However, 
this  is  generally  not  true  for  multi- 
degree-of-freedom  systems  typical  of 
finite  element  models.  The  reason 
for  this  is  that  the  vibration  modes 
associated  with  higher  frequencies 
(shorter  periods)  generally  have 
very  low  participation  factors.  Thus, 
even  though  the  higher  modes  may  be 
integrated  with  significant  error, 
their  net  contribution  is  masked  by 
the  large  participation  factors  of 
the  dominant  lower  modes  whose  in- 
tegration is  more  accurate  due  to  the 
longer  periods.  Also,  it  is  known 
that  the  higher  vibration  modes  of 
a finite  element  model  are  “ ficti- 
tious’ * in  the  sense  that  they  represent  a lumping  of  the  infinite  set  of 
higher  modes  of  vibration  corresponding  to  the  continuous  system  being 
modeled.  Therefore,  the  low  participation  of  the  fictitious  vibration 
modes  is  not  merely  a fortuitous  happenstance,  but  rather  is  a necessary 
consequence  of  proper  spatial  discretization. 

To  summarize,  for  the  cases  3 < 1 /A  the  allowable  time  step  is 
generally  governed  by  the  stability  criterion  (Equation  13  or  14)  based 
on  the  shortest  period.  Unlike  the  accuracy  consideration,  the  shortest 
period  (highest  mode)  cannot  be  ignored  because,  even  though  it  initially 
has  a low  participation  factor,  if  it  is  integrated  in  the  unstable  range, 
its  contribution  will  grow  without  bound. 

For  the  unconditionally  stable  case  3 * 1/4,  the  allowable  time  step 
is  governed  strictly  by  accuracy  considerations.  Since  the  analyst  may 
have  no  idea  of  the  frequency  content  of  the  mesh  nor  what  frequencies 
will  be  dominant  for  a specified  loading,  it  is  difficult  to  determine 
a priori  an  optimum  time  step  size.  Accordingly,  computational  experi- 
mentation is  generally  the  most  direct  method.  As  a guide,  the  time  step 
should  be  at  least  sufficiently  small  to  define  the  shape  and  character 
of  the  loading  function. 

In  the  next  section  computational  experimentation  is  performed  on 
the  test  problem  posed  in  the  objective  of  this  report. 


Figure  2.  Period  accuracy 
of  Newmark  3-method. 


COMPUTATIONAL  INVESTIGATION 

The  finite  element  code  FEAP  [9]  is  used  to  study  the  efficiency  of 
the  implicit  algorithm  3*1/4  versus  the  explicit  algorithm  3*0. 

FEAP’ s explicit  and  implicit  algorithms  accept  identical  finite  element 
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input  records  and  use  the  same  element  stiffness  routine.  Thus,  it  may 
be  presumed  the  subsequent  comparisons  provide  an  unbiased  assessment 
of  integration  efficiency. 

The  test  problem  described  is  representative  of  a class  of  problems 
of  current  interest  to  t)NA  [1,  2], 

Test  Problem  and  Finite  Element  Model 

The  test  problem  is  illustrated  in  Figure  3 which  shows  a system 
composed  of  a homogeneous  media  representative  ov.  a soft  rock-like 
material  called  tuff  and  a 3-foot  (0.914  m)  circular  steel  liner  0.78 
inch  (1.98  cm)  thick.  All  materials  are  elastic  and  are  defined  in 
Figure  3.  Plane  strain  geometry  is  assumed,  and  the  system  is  symmetric 
about  the  vertical  centerline. 

Pressure,  p(t) 


Pressure  loading  is  distributed 
uniformly  along  the  surface  of  the 
tuff,  and  the  loading  function  is  a 
triangular  pulse  whose  rise  time  is 
10  ms  with  a maximum  pressure  of  0.5 
kbar  followed  by  a 30-ms  decay  to 
zero  pressure. 

The  finite  element  mesh  topology 
representing  the  test  problem  is 
shown  in  Figure  4.  Ail  elements  are 
compatible  four-node,  isoparametric 
elements.  The  steel  liner  is  coarse- 
ly modeled  with  a single  layer  of 
32  elements  forming  a semicircle. 
Although  the  single  layer  is  insuf- 
ficient to  accurately  capture  bend- 
ing of  the  liner,  it  suffices  for 
this  investigation  where  the  concern 
is  with  time  integration, 

A simple  lumped-mass  procedure 
is  used  for  both  explicit  and 
implicit  integration  schemes. 

Boundary  conditions,  degrees  of 
freedom,  average  bandwidth,  and 
other  information  are  displayed  in 
Figure  4. 

Explicit  Results  (g  ■ 0) 

Solutions  were  attempted  for 
time  steps  of  At  • 0,0025,  0.0035, 
and  0.0045  ms.  The  solutions  for 
At  ■ 0.0025  and  0,0035  ms  were  prac- 
tically identical,  indicating  these 
solutions  provided  accurate  time 
integrations.  A typical  result  is 
shown  by  the  solid  line  in  Figure  5, 
which  is  a normalized  time-history 
plot  of  the  thrust  stress  (i.e., 
average  hoop  stress)  in  the  steel 
liner  at  the  springline. 

For  the  case  At  ■ 0,0045  ms, 
the  responses  in  the  steel  liner 
became  unstable  (wildly  erratic) 
after  approximately  150  time  steps 
or  0.68  m3.  Thereafter,  instability 
quickly  spread  throughout  the  entire 
system  in  a matter  of  a few  time 


Elements  = 381 

Degrees  freedom  = 861 
Average  bandwidth  = 25 


Loading 


Figure  4.  Finite  element 
mesh  of  test  problem. 
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steps.  Theoretically,  the  sonic 
travel  time  for  a surface  disturbance 
to  travel  through  18  feet  (5.49  m)  of 
tuff  and  excite  the  liner  is  3.8  ms. 
However,  the  numerical  wave  speed  can 
travel  vertically  downward  one  ele- 
ment dept.i  per  time  step  and  prema- 
turely excite  the  instability  of  the 
steel  liner  well  ahead  of  the  sonic 
travel  time. 

The  heuristic  stability  pre- 
diction for  maximum  allowable  At 
(Equation  14)  based  on  the  steel 
liner  elements  is  At  ■ 0.0033  ms 
[i.e.,  h ■ 0.78  inch  (1.98  cm),  c « 
235  in. /ms  (597  cm/ms)].  This  is  in 
good  agreement  with  the  observed 
instability  occurring  between  0.0035 
< At  < 0.0045  ms.  Computer  time  for 
the  central  processor  is  summarized 
in  Table  1. 

Implicit  Results  (8  ■ 1/4) 


Solutions  were  obtained  for  the  Flgu“  ?-,"1St°5y  plots  of 
time  steps  it  - 0.4,  0.8,  and  2.0  ms,  liner  thrust  stress, 

uhlch  are  two  orders  of  magnitude  Implicit  and  explicit, 

p' eater  than  the  allowable  explicit 

time  steps.  Since  B * 1/4  provides  an  unconditionally  stable  implicit 
algorithm,  the  only  concern  is  with  accuracy.  In  each  case,  the  results 
were  in  close  agreement  with  each  other  and  the  explicit  solution  as 
illustrated  by  the  discrete  symbols  superimposed  on  the  explicit  solution 
presented  in  Figure  5.  A slight  error  can  be  observed  near  the  peak  of 
the  response  for  the  case  At  * 2.0  ms.  Accordingly,  for  the  sake  of 
subsequent  comparisons,  it  will  be  said  that  the  maximum  allowable  time 
step  based  on  accuracy  is  At  ■ 0.8  ms.  This  corresponds  to  12  time  steps 
within  the  loading  rise  time.  Computer  cost  (time)  for  the  central 
processor  is  summarized  in  Table  1. 


Efficiency  of  Explicit  Versus  Implicit 


The  computer  times  reported  in  Table  1 were  determined  by  assessing 
a system  “clock  routine,* * which  measures  the  actual  running  time  of 
each  step  for  both  the  explicit  and  implicit  algorithms.  These  results 
indicate  the  number  of  time  steps  required  to  complete  the  40  ms  loading 
duration  is  11500  (i.e,,  40/0.0035)  for  the  explicit  case,  but  only  50 
(i.e.,  40/0.8)  for  the  implicit  case;  or  a ratio  of  230  in  favor  of  the 
implicit  method.  On  the  other  hand,  the  computer  cost  (time)  per  step  of 
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Table  1. 


Central  Processor  Time  for  Explicit 
and  Implicit  Algorithms 


Computer  Time 

Computer  Time 

Algorithm 

Per  Step  of  Algorithm 

Per  Complete  Solution 

(s) 

(s) 

Explicit 
(8  » 0) 

0.16A 

1,870.0 

Implicit 
(6  - 1/A) 

A.  2325  or  1.22C 

6A.0 

i 

aComplete  solution  time  is  based  on  AO  ms  problem,  using 
At  « 0.0035  ms  for  explicit  and  At  ■ 0.8  ms  for  implicit. 

^Time  for  a solution  step  requiring  triangularization. 

cTime  for  a solution  step  not  requiring  triangularization. 


the  implicit  method,  where  the  stiffness  matrix  is  triangularized,  is  25 
times  more  than  an  explicit  step  cost.  For  subsequent  steps  not  requiring 
triangularization,  the  implicit  cost  is  8 times  more  than  the  explicit 
cost  per  step.  As  a net  result  the  computer  cost  (based  on  central 
processing  time)  of  the  explicit  solution  is  30  times  more  expensive 
than  the  implicit  solution  for  a complete  AO-ms  run. 

The  observed  efficiency  of  the  implicit  method  is  primarily  due  to 
the  linear  nature  of  the  test  problem. 

Nonlinear  systems  do  not  appreciably  alter  the  computer  costs  for 
explicit  algorithms  because  stiffness  changes  are  dealt  with  at  the 
element  level  on  the  right-hand  side  of  the  equilibrium  equations. 
Moreover,  since  the  explicit  step  size  is  inherently  small,  iterations 
within  the  time  step  are  generally  not  required.  Conversely,  nonlinear 
stiffnesses  in  the  implicit  method  may  require  triangularizing  and 
iterating  within  every  time  step,  as  well  as  reducing  the  size  of  the 
time  step.  For  example,  previous  experience  indicates  that,  if  the  tuff 
material  is  modeled  with  a nonhardening  plasticity  law,  tl>ui  to  maintain 
accuracy  for  the  implicit  algorithm  requires  At  = 0.2  ms,  as  well  as 
triangularizing  and  iterating  within  each  time  step.  Under  these  condi- 
tions, the  computer  cost  can  be  estimated  from  Table  1 to  be  about 
equivalent  to  the  cost  of  the  explicit  method. 

Summarizing  for  the  class  of  problems  considered,  the  implicit 
algorithm  is  significantly  more  efficient  for  linear  systems;  however, 
for  nonlinear  systems  the  two  methods  are  competitive.  The  most  prom- 
ising efficient  scheme  is  a combined  explicit/implicit  algorithm 
(presented  next). 
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COMBINING  EXPLICIT  AND  IMPLICIT 

The  nature  of  the  foregoing  boundary  value  problem  provides  an 
illustrative  example  of  the  potential  advantages  of  combining  the 
explicit  and  implicit  methods  into  a single  algorithm.  To  see  this, 
consider  the  advantages  of  applying  the  implicit  algorithm  to  nodes 
(degrees  of  freedom)  associated  with  the  steel  liner,  while  the  remain- 
ing nodes  (degrees  of  freedom)  in  the  tuff  are  integrated  with  an 
explicit  algorithm.  With  this  assignment,  explicit  stability  is 
governed  by  the  larger  and  sonically  slower  tuff  elements,  resulting  in 
a stable  time  step  on  the  order  of  10  to  100  times  larger  than  the  time 
step  required  when  the  steel  elements  govern.  By  the  same  token,  the 
implicit  algorithm  is  significantly  enhanced  because  only  the  small 
portion  of  the  stiffness  matrix  associated  with  the  steel  liner  needs 
triangularization. 

The  above  concepts  are  presented  more  formally  in  the  following 
general  discussion.  Consider  an  arbitrary  body  discretized  by  finite 
elements  such  that  all  * *i**  nodes  are  to  be  integrated  implicitly  and 
all  “e”  nodes  are  to  be  integrated  explicitly.  These  nodes  are  denoted 
in  Figure  6 by  a dashed  line  encompassing  the  “i”  nodes,  and  all  nodes 
exterior  to  the  dashed  curve  are  “e?’  nodes.  Those  element  stiffnesses 
associated  with  only  ‘<1**  nodes  are  denoted  by  Kj.i,  while  those  asso- 
ciated with  only  “e”  nodes  are  expressed  as  Eeg.  Mixed  elements 
associated  with  both  “e1  * and  “i» » nodes  are~denoted  by  gei. 

With  these  definitions,  the  equilibrium  equation  (Equation  1)  can 
be  written  in  partitioned  form  as: 


Kee  I Kei 
*i"e‘  J hi 

• I • 


(17) 


where 


K 


ee 


K 


ii 


K 


ei 


displacements  at  explicit,  implicit  nodes 

accelerations  at  explicit,  implicit  nodes 

lumped  masses  at  explicit,  implicit  nodes 

nodal  forces  at  explicit,  implicit  nodes 

global  stiffness  of  elements  with  only  **e* ’ nodes 

global  stiffness  of  elements  with  only  “i’»  nodes 
T 

K^e  « global  stiffness  of  elements  with  both  e and  i 
nodes 
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Figure  6.  Conceptualization  of  explicit  and  implicit 
integration  regions. 


Expressing  the  partitioned  equilibrium  equation  as  two  coupled  systems 
of  equations  gives 


[Me]  {u,} 

+ tKee] 

tue) 

+ 'Kei] 

{u^}  - 

<£e> 

(18) 

(MJ  {\1±} 

+ [Ku] 

+ (Kle! 

tue)  - 

(19) 

To  numerically  integrate  Equations  18  and  19  step-by-step,  it  is  assumed 
the  responses  {ug}  , {u^} and  their  first  two  derivatives  are  known  for 
time  t,  and  the  objective  is  to  determine  the  responses  {u  )t+At»^ui^t+At 
and  their  first  two  derivatives  for  time  t + At.  To  this  end,  the  follow- 
ing procedure  constitutes  a combined  explicit-implicit  algorithm  using 
the  Newmark  B-method  defined  by  Equations  2 and  3. 

Step  1.  Predict  the  **e* * nodal  displacements  directly 

from  Equation  3 with  8 » 0,  (i.e,,  explicit)  to  get: 


{u.W 
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where 


Step  2.  Again  using  Equation  3 with  B 'f  0,  express  fu^t^t  as  a 
function  of  {«•(.’ t+At  and  previous  responses  to  get: 
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where  {a.^ 
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Step  3.  Inserting  the  above  expression  for  {ui)t+At  into  Equation 
19,  (ui}t+At  is  determined  by  solving  the  coupled  system: 
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Step  4.  Having  determined  and  {ui}t+At>  compute  the  ^ 

* corresponding  accelerations  using  Equation  18  for 

and  the  first  equation  in  Step  2 for  (u^t+At  t0  Set: 
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Step  5.  Lastly,  update  the  velocities  using  Newroark’s  expression, 
Equation  2 
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This  procedure  is  repeated  step-by-step  throughout  the  time  of 
interest. 

The  advantages  of  the  combined  explicit-implicit  algorithm  is  quite 
remarkable  in  that  it  is  potentially  one  or  two  orders  of  magnitude 
more  efficient  than  either  method  applied  individually.  Note  that  the 
only  triangularization  required  is  in  Step  3 for  the  submatrix  [%*] 
whose  size  is  small  compared  to  the  entire  system  (e.g.,  steel  elements). 
Moreover,  the  average  bandwidth  of  [K^i]  is  easily  minimized  by  judicious 
compact  node  numbering  of  the  “i”  nodes,  which  implies  rapid  solutions 
with  minimal  storage.  Node  numbering  of  ‘*e*’  nodes  has  no  bearing  on 
the  solution  efficiency  since  all  stiffness  operations  involving  “e** 
nodes  are  matrix  multiplications  perform,  •!  at  the  element  level.  Since 
the  domain  of  the  “e**  nodes  are  selected  to  permit  the  explicit 
stability  criterion  to  be  based  -on  the  larger  and  sonically  slower 
elements,  the  critical  time  step  may  be  increased  by  one  or  two  orders 
of  magnitude,  depending  on  the  nature  of  the  problem. 

It  is  easy  to  conceive  of  large,  three-dimensional,  nonlinear 
problems  such  that  a combined  method  as  presented  here  may  be  the  only 
economically  feasible  alternative  of  obtaining  a solution.  Hence,  com- 
bined methods  should  be  vigorously  investigated. 


CONCLUSIONS  AND  RECOMMENDATIONS 

Within  the  class  of  problems  and  limitations  defined  in  this  study, 
the  following  comments  appear  valid  for  numerical  integration  by  the 
Newmark  6-method. 

1.  The  maximum  allowable  time  step  for  conditionally  stable  schemes 

(6  < 1/4)  is  generally  controlled  by  the  stability  criterion,  Equation  13 
or  14, 

2.  The  maximum  allowable  time  step  for  the  unconditionally  stable 
scheme  (6  ■ 1/4)  is  controlled  by  accuracy  considerations  which  may  be 
determined  by  computational  experimentation.  As  a guide,  define  the 
time  step  to  subdivide  the  rise  time  into  10  increments. 

3.  For  the  linear  system  studied  and  computer/program  employed,  the 
implicit  algorithm  (6  * 1/4)  allows  a time  step  more  than  200  times 
larger  than  the  explicit  method  (6  * 0)  with  equivalent  accuracy.  How- 
ever, each  step  of  the  implicit  algorithm  is  8 to  25  times  more  costly 
in  computer  time  than  a step  of  the  explicit  algorithm  depending  on 
whether  or  not  triangularization  is  required.  As  a net  result  the 
implicit  method  was  approximately  30  times  more  efficient. 

4.  Nonlinear  systems  penalize  implicit  algorithms  much  more  severely 
than  explicit  algorithms  fo  the  extent  that  the  two  methods  become 
equally  competitive  for  the  class  of  problems  considered. 
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5.  The  combined  explicit/implicit  algorithm  has  the  potential  of 
enhancing  efficiency  by  one  or  more  orders  of  magnitude.  It  is  recom- 
mended that  combined  integration  methods  be  thoroughly  investigated. 
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CORNELL  UNIVERSITY  Ithaca  NY  (Serials  Dept.  Engr  Lib.) 

DAMES  ti  MOORE  LIBRARY  LOS  ANGELES,  CA 

DUKE  UNIV  MEDICAL  CENTER  B.  Mtiga,  Durham  NC;  DURHAM,  NC(VESIC) 

FLORIDA  ATLANTIC  UNIVERSITY  BOCA  RATON.  FL (MC  ALLISTER);  Boca  Raton  FL  (Ocean  Engr  Dept.,  C, 
Lin) 

FLORIDA  ATLANTIC  UNIVERSITY  Boca  Raton  FL  (W.  Tessin) 

FLORIDA  TECHNOLOGICAL  UNIVERSITY  OR  LAN  IX),  FL  (HARTMAN) 

GEORGIA  INSTITUTE  OF  TECHNOLOGY  Atlanta  GA  (School'of  Civil  Engr.,  Kahn);  Atlanta  GA  (B.  Mazanli) 
IOWA  STATE  UNIVERSI  TY  Ames  IA  (CE  Dept,  Handy) 

VIRGINIA  INST.  OF  MARINE  SCI.  Gloucester  Point  VA  (Library) 

LEHIGH  UNIVERSITY  BETHLEHEM.  PA  (MARINE  GEOTECHNICAL  LAB.,  RICHARDS):  Bethlehem  PA 
(Fritz,  Engr.  Lab  No.  13,  Beedlc);  Bethlehem  PA  (Lindcrman  Lib.  No.30,  Flcckstcincr) 

LIBRARY  OF  CONGRESS  WASHINGTON.  DC  (SCIENCES  Si  TECH  DIV) 

MASSACHUSETTS  INST.  OF  TECHNOLOGY  Cambridge  MA  (Rm  10-500,  Tech.  Reports.  Engr.  Lib.);  Cambridge 
MA  (Rm  14  E210,  Tech.  Report  Lib.);  Cambridge  MA  (Whitman) 

MICHIGAN  TECHNOLOGICAL  UNIVERSITY  Houghton.  Ml  (Haas) 

NORTHWESTERN  UNIV  Z.P.  Bazant  Evanston  IL 

NY  CITY  COMMUNITY  COLLEGE  BROOKLYN.  NY  (LIBRARY) 

OREGON  STATE  UNIVERSI  TY  CORVALLIS,  OR  (CE  DEFT,  BELL);  CORVALLIS,  OR  (CE  DEPT,  HICKS) 
PENNSYLVANIA  STATE  UNIVERSI  TY  UNIVERSITY  PARK,  PA  (GOTOLSKI) 

PURDUE  UNIVERSITY  Lafayette  IN  (Leonards);  Lafayette,  IN  (Allsehaeffl);  Lafayette,  IN  (CE  Engr.  Lib) 

SAN  DIEGO  STATE  UNIV.  Dr.  Krishnamoorlhy,  San  Diego  CA 
SOUTHWEST  RSCH  INST  J.  Matson,  San  Antonio  TX;  R.  DeHart,  San  Antonio  TX 
STANFORD  UNIVERSITY  Stanford  CA  (Gene) 

STATE  UNIV.  OF  NEW  YORK  Buffalo,  NY 

TEXAS  A&M  UNIVERSITY  COLLEGE  STATION,  TX  (CE  DEPT);  College  TX  (CE  Dept.  Herbich) 
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U.S.  MERCHANT  MARINE  ACADEMY  KINGS  POINT.  NY  (REPRINT  CUSTODIAN) 

UNIVERSITY  OF  CALIFORNIA  BERKELEY.  CA  (CE  DEPT,  GERWICK);  BERKELEY.  CA  (CE  DEPT. 

MITCHELL):  BERKELEY.  CA  (OFF.  BUS.  AND  FINANCE.  SAUNDERS):  Berkeley  CA  (B.  Brcslcr);  DAVIS, 
CA  (CE  DEPT.  TAYLOR);  LIVERMORE.  CA  (LAWRENCE  LIVERMORE  LAB,  TOKARZ);  La  Jolla  CA  (Acq. 
Dept.  Lib.  C-075A):  Los  Angeles  CA  (Engr  I.  K.  Lee);  SAN  DIEGO,  CA,  LA  JOLLA,  CA  (SEROCKI) 
UNIVERSITY  OF  DELAWARE  Newark,  DE  (Dept  of  Civil  Engineering,  Chesson) 

UNIVERSITY  Of  HAWAII  HONOLULU.  HI  (CE  DEPT,  GRACE):  HONOLULU.  HI  (SCIENCE  AND  TECH. 
DIV.);  Honolulu  HI  (Dr.  Szilard) 

UNIVERSITY  OF  ILLINOIS  Metz  Ref  Rrn.  Urbana  IL;  URBANA,  IL  (DAVISSON);  URBANA,  IL  (LIBRARY), 
URBANA,  IL  (NEWARK):  Urbana  IL  (CE  Dept,  W.  Gamble) 

UNIVERSITY  OF  MASSACHUSETTS  (Heroncmus).  Amherst  MA  CE  Dept 
UNIVERSITY  OF  MICHIGAN  Ann  Arbor  Ml  (Richart) 

UNIVERSITY  OF  NEBRASKA-LINCOLN  Lincoln.  NE  (Ross  Ice  Shelf  Proj.) 

UNIVERSITY  OF  NEW  MEXICO  Albuquerque  NM  (Soil  Mcch.  & Pav.  Div..  J.  Nielsen) 

UNIVERSITY  OF  TEXAS  Inst.  Marine  Sci  (Library).  Port  Aransas  TX 

UNIVERSITY  OF  TEXAS  AT  AUSTIN  AUSTIN.  TX  (THOMPSON):  Austin  TX  (R.  Olson):  Austin,  TX  (Breen) 
UNIVERSITY  OF  WASHINGTON  Dept  of  Civil  Engr  (Dr.  Mattock).  Seattle  WA:  SEATTLE.  WA  (MERCHANT); 

SEATTLE.  WA  (OCEAN  ENG  RSCH  LAB.  GRAY);  Seattle  WA  (E.  Linger) 

URS  RESEARCH  CO.  LIBRARY  SAN  MATEO.  CA 

US  GEOLOGICAL  SURVEY  Off.  Marine  Geology.  Mailstop  915.  Reston  VA 

AEROSPACE  CORP.  Acquisition  Group,  Los  Angeles  CA 

ALFRED  A.  YEE  & ASSOC.  Honolulu  HI 

ARVID  GRANT  OLYMPIA,  WA 

ATLANTIC  RICHFIELD  CO.  DALLAS. TX  (SMITH) 

AUSTRALIA  Dept.  PW  (A.  Hicks),  Melbourne 
BECHTEL  CORP.  SAN  FRANCISCO.  CA  (PHELPS) 

BELGIUM  NAECON,  N.V.,  GEN. 

BET  HLEHEM  STEEL  CO.  BETHLEHEM,  PA  (ST  EELE) 

BROWN  & ROOT'  Houston  TX  (D.  Ward) 

CANADA  Mein  Univ  Newfoundland  (Chari),  St  Johns;  Surveyor,  Nenningcr&  Chencvert  Inc..  Montreal;  Warnoek 
Hcrsey  Prof.  Srv  Ltd,  La  Sale,  Quebec 
CF  BRAUN  CO  Du  Bouchct,  Murray  Hill,  NJ 
CONCRETE  TECHNOLOGY  CORP.  TACOMA,  WA  (ANDERSON) 

CONRAD  ASSOC.  Van  Nuys  CA  (A.  Luisoni) 

DRAVO  CORP  Pittsburgh  PA  (Giannino);  Pittsburgh  PA  (Wright) 

NORWAY  DET  NORSKE  VERITAS  (Library).  Oslo 
EVALUATION  ASSOC.  INC  KING  OF  PRUSSIA.  PA  (FEDELE) 

FORD,  BACON  & DAVIS.  INC.  New  York  (Library) 

FRANCE  Dr.  Dutcrlrc,  Boulogne;  P.  Jensen,  Boulogne 
GEOTECHNICAL  ENGINEERS  INC.  Winchester,  MA  (Paulding) 

GL1DDEN  CO.  STRONGSVILLE.  OH  (RSCH  LIB) 

GLOBAL  MARINE  DEVELOPMENT  NEWPORT  BEACH,  CA  (HOLLETT) 

GRUMMAN  AEROSPACE  CORP.  Bcthpage  NY  (Tech.  Info.  Ctr) 

HALEY  & ALDRICH,  INC.  Cambridge  MA  (Aldrich.  Jr.) 

HONEYWELL,  INC.  Minneapolis  MN  (Residential  Engr  Lib.) 

HUGHES  AIRCRAFT  Culver  City  CA  (Tech.  Doc.  Ctr) 

ITALY  M.  Caironi,  Milan;  Sergio Tattoni  Milano 
JAMES  CO.  R.  Girdley,  Orlando  FL 

LAMONT-DOHERTY  GEOLOGICAL  OBSERV.  Palisades  NY  (McCoy);  Palisades  NY  (Sebvyn) 

LOCKHEED  MISSILES  & SPACE  CO.  INC.  SUNNYVALE,  CA  (PHILLIPS);  Sunnyvale  CA  (Rynewicz) 
LOCKHEED  OCEAN  LABORATORY  San  Diego  CA  (F.  Simpson) 

MARATHON  OIL  CO  Houston  TX  (C.  Seay) 

MCCLELLAND  ENGINEERS  INC  Houston  TX  (B.  McClelland) 

MCDONNEL  AIRCRAFT  CO.  Dept  501  (R.H.  Fayman),  St  Louis  MO 
MEDALL  & ASSOC.  INC.  J.T.  GAFFEY  II  SANTA  ANA.  CA 
MEXICO  R.  Cardenas 

MOBIL  PIPE  LINE  CO.  DALLAS,  TX  MGR  OF  ENGR  (NOACK) 

MUESER.  RUTLEDGE.  WENT  WORTH  AND  JOHNSTON  NEW  YORK  (RICHARDS) 

NEW  ZEALAND  New  Zealand  Concrete  Research  Assoc.  (Librarian),  Porirua 
NEWPORT  NEWS  SH1PBLDG  & DRYDOCK  CO.  Newport  News  VA  (Tech.  Lib.) 
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NORWAY  DET  NORSKE  VERITAS  (Roren)  Oslo:  J.  Creed.  Ski:  Norwegian  Tech  Univ  (Brandtzacg),  Trondheim 
OEESHORE  DEVELOPMENT  ENG.  INC.  BERKELEY.  CA 
PACIFIC  MARINE  TECHNOLOGY  LONG  BEACH.  CA  (WAGNER) 

PORTLAND  CEMENT  ASSOC.  SKOKIE.  IL (CORELY):  Skokie  IL(Rsch  & Dev  Lab,  Lib.) 

PRESCON  CORPTOWSON,  MD  (KELLER) 

RAND  CORP.  Santa  Monica  CA  (A.  Lanpa) 

RAYMOND  INTERNATIONAL  INC.  CHERRY  HILL.  NJ  (SOILTECH  DEPP) 

RIVERSIDE  CEMENT  CO  Riverside  CA  (W.  Smith) 

SANDIA  LABORATORIES  Library  Div..  Livermore  CA 
SCHUPACK  ASSOC  SO.  NORWALK.  CT(SCHUPACK) 

SEATECH  CORP.  MIAMI.  FL  (PERONI) 

SHELL  DEVELOPMENT  CO.  Houston TX  (E.  Doyle) 

SHELL  OIL  CO.  HOUSTON,  TX  (MARSHALL) 

SWEDEN  GcoTcch  Inst:  VBB  (Library),  Stockholm 
TIDEWATER  CONSTR.  CO  Norfolk  VA  (Fowler) 

TRW  SYSTEMS  CLEVELAND.  OH  (ENG.  LIB.);  REDONDO  BEACH.  CA  (DAI) 

UNITED  KINGDOM  Cement  & Concrete  Assn.  (R  Rowe)  Wcxhatn  Springs.  Slough  Buck:  Cement  & Concrete 
Assoc.  (Lit,  Ex).  Bucks;  D.  New,  G.  Maunsclt  & Partners,  London;  Shaw  & Hatton  (F.  Hansen),  London;  Taylor, 
Woodrow  Constr  (0I4P),  Southall,  Middlesex;  Univ.  of  Bristol  (R.  Morgan).  Bristol 
USGS  MENLO  PARK.  CA  (YOUD) 

WESTINGHOUSE  ELECTRIC  CORP,  Annapolis  MD  (Oceanic  Div  Lib,  Bryan);  Library,  Pittsburgh  PA 
WISS,  JANNEY,  ELSTNER.  & ASSOC  Northbrook.  IL  (J.  Hanson) 

WOODWARD-CLYDE  CONSULTANTS  Oakland  CA  (A.  Harrigan):  PLYMOUTH  MEETING  PA  (CROSS,  III) 

AL SMOOTS  Los  Angeles.  CA 
BULLOCK  La  Canada 
F.  HEUZE  Boulder  CO 
HAMEED  ELNAGGAR  Wexford  PA 
CAPT  MURPHY  SAN  BRUNO,  CA 
GREG  PAGE  EUGENE.  OR 
R.F.  BUSIER  Old  Saybrook  CT 
T.W,  MERMEL  Washington  DC 
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